Learning to open HDF5 files in R

# load libraries
library(raster)
## Loading required package: sp
library(rhdf5)
library(rgdal)
## rgdal: version: 1.1-10, (SVN revision 622)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-3

Let’s open a file

# set file path
f <- "../NEONdata/D17-California/TEAK/2013/spectrometer/reflectance/Subset3NIS1_20130614_100459_atmcor.h5"

# view h5 structure
h5ls(f)
##    group                                 name       otype  dclass
## 0      /                     ATCOR_Input_File H5I_DATASET  STRING
## 1      /                 ATCOR_Processing_Log H5I_DATASET  STRING
## 2      /                Aerosol Optical Depth H5I_DATASET INTEGER
## 3      /                               Aspect H5I_DATASET   FLOAT
## 4      /                          Cast Shadow H5I_DATASET INTEGER
## 5      / Dark Dense Vegetation Classification H5I_DATASET INTEGER
## 6      /                 Haze-Cloud-Water Map H5I_DATASET INTEGER
## 7      /                  Illumination Factor H5I_DATASET INTEGER
## 8      /                          Path Length H5I_DATASET   FLOAT
## 9      /                   Processing Version H5I_DATASET  STRING
## 10     /                          Reflectance H5I_DATASET INTEGER
## 11     /                Shadow_Processing_Log H5I_DATASET  STRING
## 12     /                      Sky View Factor H5I_DATASET INTEGER
## 13     /               Skyview_Processing_Log H5I_DATASET  STRING
## 14     /                                Slope H5I_DATASET   FLOAT
## 15     /          Slope_Aspect_Processing_Log H5I_DATASET  STRING
## 16     /                  Solar Azimuth Angle H5I_DATASET   FLOAT
## 17     /                   Solar Zenith Angle H5I_DATASET   FLOAT
## 18     /                    Surface Elevation H5I_DATASET   FLOAT
## 19     /                 Visibility Index Map H5I_DATASET INTEGER
## 20     /                   Water Vapor Column H5I_DATASET INTEGER
## 21     /             coordinate system string H5I_DATASET  STRING
## 22     /                       flightAltitude H5I_DATASET   FLOAT
## 23     /                        flightHeading H5I_DATASET   FLOAT
## 24     /                           flightTime H5I_DATASET   FLOAT
## 25     /                                 fwhm H5I_DATASET   FLOAT
## 26     /                             map info H5I_DATASET  STRING
## 27     /              to-sensor azimuth angle H5I_DATASET   FLOAT
## 28     /               to-sensor zenith angle H5I_DATASET   FLOAT
## 29     /                           wavelength H5I_DATASET   FLOAT
##                dim
## 0                1
## 1                1
## 2    544 x 578 x 1
## 3    544 x 578 x 1
## 4    544 x 578 x 1
## 5    544 x 578 x 1
## 6    544 x 578 x 1
## 7    544 x 578 x 1
## 8    544 x 578 x 1
## 9                1
## 10 544 x 578 x 426
## 11               1
## 12   544 x 578 x 1
## 13               1
## 14   544 x 578 x 1
## 15               1
## 16           1 x 1
## 17           1 x 1
## 18   544 x 578 x 1
## 19   544 x 578 x 1
## 20   544 x 578 x 1
## 21               1
## 22         5732053
## 23         5732053
## 24         5732053
## 25         426 x 1
## 26               1
## 27   544 x 578 x 1
## 28   544 x 578 x 1
## 29         426 x 1

Import spatial information

# import spatial info
mapInfo <- h5read(f, 
                  "map info",
                  read.attributes = TRUE)
mapInfo
## [1] "UTM,1,1,325963.0,4103482.0,1.0000000000e+000,1.0000000000e+000,11,North,WGS-84,units=Meters"
## attr(,"Description")
## [1] "Basic Map information for envi style programs"

Grab reflectance metadata

# read in reflectance data attributes
reflInfo <- h5readAttributes(f, "Reflectance")
reflInfo
## $DIMENSION_LABELS
## [1] "Wavelength" "Line"       "Sample"    
## 
## $Description
## [1] "Atmospherically corrected reflectance."
## 
## $`Scale Factor`
## [1] 10000
## 
## $Unit
## [1] "unitless. Valid range 0-1."
## 
## $`data ignore value`
## [1] "15000.0"
# define scale factor
scaleFactor <- reflInfo$`Scale Factor`

# define no data value
noDataValue <- reflInfo$`data ignore value`  # it's a character
noDataValue <- as.numeric(noDataValue)
str(noDataValue)  # looks good
##  num 15000

Import data dimensions

# open file for viewing
fid <- H5Fopen(f)

# open the reflectance dataset
did <- H5Dopen(fid, "Reflectance")
did  # note dimensions are in columns x rows x bands
## HDF5 DATASET
##         name /Reflectance
##     filename 
##         type H5T_STD_I16LE
##         rank 3
##         size 544 x 578 x 426
##      maxsize 544 x 578 x 426
# grab dataset dimensions
sid <- H5Dget_space(did)
dims <- H5Sget_simple_extent_dims(sid)$size
dims  # columns, rows, bands; r reads rows 1st, columns 2nd
## [1] 544 578 426
# close all open connections
H5Sclose(sid)
H5Dclose(did)
H5Fclose(fid)

Read in reflectance data

# extract slice of H5 file
b56 <- h5read(f,
              "Reflectance",
              index = list(1:dims[1], 1:dims[2], 56))
b56

class(b56)

Convert data to matrix

# convert to matrix
b56 <- b56[,,1]  # z dimension goes away, set to 1
class(b56)
## [1] "matrix"
# let's plot some data
image(b56)

# stretch the image by log transforming
image(log(b56), main = "Log transformed data")

# force non-scientific notation
options("scipen"=100, "digits"=4)

# look at histograms
hist(b56)

# histogram of stretched image
hist(log(b56))

Time to clean up our data

Note that data is stored as integers rather than floats. This way, data takes up less space. Can use scale factor to convert to decimal point at later time.

# assign no data values to object
b56[b56 == noDataValue] <- NA

# apply scale factor
b56 <- b56 / scaleFactor  # grabbed scale factor earlier
hist(b56)  # right skew

Transpose data

# transpose data, flip rows and columns
b56 <- t(b56)
image(log(b56))  # still not quite right

Create spatial extent

# split out map info object
mapInfo <- strsplit(mapInfo, ",")
mapInfo <- unlist(mapInfo)  # remove nesting
mapInfo
##  [1] "UTM"               "1"                 "1"                
##  [4] "325963.0"          "4103482.0"         "1.0000000000e+000"
##  [7] "1.0000000000e+000" "11"                "North"            
## [10] "WGS-84"            "units=Meters"
# value at 3 element in list
mapInfo[3]  # note that it's a character
## [1] "1"
# define upper left hand corner coordinate
xMin <- as.numeric(mapInfo[4])
yMax <- as.numeric(mapInfo[5])

# get spatial resolution
xRes <- as.numeric(mapInfo[6])
yRes <- as.numeric(mapInfo[7])

# define lower right hand corner coordinate
xMax <- xMin + (dims[1] * xRes)
yMin <- yMax - (dims[2] * yRes)

# create extent object
rasExt <- extent(xMin, xMax, yMin, yMax)
rasExt
## class       : Extent 
## xmin        : 325963 
## xmax        : 326507 
## ymin        : 4102904 
## ymax        : 4103482
# create raster object
b56r <- raster(b56, 
               crs = CRS("+init=epsg:32611"))

extent(b56r) <- rasExt
b56r
## class       : RasterLayer 
## dimensions  : 578, 544, 314432  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : 325963, 326507, 4102904, 4103482  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:32611 +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : 0.0033, 0.5678  (min, max)
# plot data
plot(b56r, main="Spatially referenced data")

Import NEON functions

# install devtools
# install.packages("devtools")
library(devtools)

# install_github("lwasser/neon-aop-package/neonAOP")
library(neonAOP)

b55 <- open_band(f,
                 bandNum = 55,
                 epsg = 32611)
plot(b55)

# import several bands
bands <- c(58, 34, 19)  # decreasing numeric order intentional to make r, g, b
bands
## [1] 58 34 19
# create raster stack
RGBStack <- create_stack(f,
                         bands = bands,
                         epsg = 32611)
plot(RGBStack)  # plot 3 bands separately

# plot RGB image
plotRGB(RGBStack,
        stretch = "lin")  # make sure ordered correctly

# cir image
bands <- c(90, 34, 19)

# create raster stack
CIRStack <- create_stack(f,
                         bands = bands,
                         epsg = 32611)

# plot RGB image
plotRGB(CIRStack,
        stretch = "lin")  # make sure ordered correctly